% Declare the relevant parameters
h=0.025;
w0=.0;
tmax=1.0;

nmax=tmax/h;

wlast=w0;
tlast=0;

t(1) = 0.0;
w(1) = 0.0;
y(1) = 0.0;

% Beginning of Midpoint Integration

for i=2:(nmax+1)
%
% Use call to fuction deriv instead of following 2 Lines !!
% w(i) = wlast + h**2*(1/2*exp(3*(tlast+h/2))+tlast*exp(3*tlast)-2*wlast)
% w(i) = w(i) + h*(tlast*exp(3*(tlast+h/2))-2wlast)
%
  a1 = wlast + h/2*deriv(tlast,wlast);
  w(i) = wlast + h*deriv(tlast + h/2,a1);
  t(i) = tlast + h;
  y(i) = exact(t(i));
  d(i) = abs(y(i) - w(i));
  wlast=w(i);
  tlast=t(i);
end 

% Plot Results for w, y

plot(t,w,'ob');
hold on;
plot(t,y,'r');

% Label the plot 

title('Solution to dy/dt=t*exp(3*t)-2y');
xlabel('Time');
ylabel('y');
legend('numerical','exact')

% Save the results (for each timestep) in a file

file = zeros(size(t,2),4);

file(:,1)=t';
file(:,2)=w';
file(:,3)=y';
file(:,4)=d';

save output file -ASCII;
