%  Euler's method to solve system of IVP's
%      dx/dt = f(t,x)
%      x(0) = x_0
%  from t=0 to t=tend.
%
%  Where x = [x_1, x_2, ..., x_n]^T
%
%  f(t,x) is specified through frhsv.m

clear;
ndim = input('enter the number of equations: ');
tend=input('enter the value of tend: ');
n=input('enter the number of steps: ');

dt=tend/n;

t = 0:dt:tend;
x = zeros(n+1,ndim);   % approximate solution from Euler's

% Read the initial conditions:
for j=1:ndim
    x(1,j) = input('enter initial values: ');
end

% Main loop:
for k=1:n
    x(k+1,:) = x(k,:) + dt*frhsv(t(k),x(k,1:ndim));
end

% The following printout assumes ndim=2.

fprintf('\n*** Euler method printout ***\n\n');

for k=1:n+1
fprintf('t = %6.4f  x_1 = %8.3f  x_2 = %8.3f\n',t(k),x(k,1),x(k,2));
end

plot(t,x(:,1),'-b',t,x(:,2),'r');
xlabel('time'); ylabel('x'); title('solution to the system IVP');
%axis([xbegin xend -6 10]);

legend('x_1','x_2');
