% Example of using ode45

t0 = 0;
tf = 15;
y0 = [1 2 1];

[t, y] = ode45('odefun', [t0 tf], y0);

plot(t, y(:, 1), 'b-', t, y(:, 2), 'g--', t, y(:, 3), 'r:')
pause

% Compare with the exact solution (known in this case)
tt = linspace(t0, tf, 100);
y1 = cos(tt) + sin(tt) - exp(-tt) + 1;
y2 = cos(tt) - sin(tt) + exp(-tt);
y3 = cos(tt) - sin(tt);
hold on
plot(tt, y1, 'bo', tt, y2, 'gx', tt, y3, 'r.')
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% function yp = odefun(t, y)
% example of a system of 3 first order ODEs
% y1' = y2;
% y2' = -y1 - 2*exp(-t) + 1;
% y3' = -y1 -   exp(-t) + 1;
function yp = odefun(t, y)

yp(1, 1) =  y(2);
yp(2, 1) = -y(1) - 2*exp(-t) + 1;
yp(3, 1) = -y(1) -   exp(-t) + 1;
