% the program takes the intial data x_0 and y_0 at the initial time t_0 and 
% computes up to time t_final using timesteps of size dt.  It creates
% three vectors, t, x and y, which you can then plot against each other
% plot(x,y), hold on, plot(x(1),y(1),'o')
%
% [t,x,y] = sys_taylor_2(x_0,y_0,t_0,t_final,dt)

function [t,x,y] = sys_taylor_2(x_0,y_0,t_0,t_final,dt)

% number of steps to take:
n = ceil((t_final-t_0)/dt) + 1;

% put initial data into the vectors
t(1) = t_0;
x(1) = x_0;
y(1) = y_0;
for i=1:n-1

% define the quantities you need to do the time-stepping for the system
%     dx/dt = x + 4 y
%     dy/dt = x + y
  xp = x(i) + 4*y(i);
  yp = x(i) + y(i);
  xpp = xp + 4*yp;	
  ypp = xp + yp;

% take the time-steps 	
  x(i+1) = x(i) + dt*xp + dt^2/2*xpp;
  y(i+1) = y(i) + dt*yp + dt^2/2*ypp;
  t(i+1) = t_0 + i*dt;
end