%Huen's method - second-order Runge Kutta
%Sove ODE y'=-2*t*y^2, time step and duration of interation
%are input paramethers.
% plot solution y(t)

function [] = huen (h,time)
y0=1;	% initial condition
xp= linspace(0,time,100); %always take 100 points
yp= y0./(1+xp.^2); % evaluate exact solution
figure(1);
plot(xp,yp,'b.')
title('exact y(t) ---   approximate y(t) ...')

yn(1)=y0;
xn= 0:h:time;
for j=1:1:time/h   % use the requested stepsize h
k1= -2*xn(j)*yn(j)^2;
k2= -2*xn(j+1)*(yn(j)+ h*k1)^2;
yn(j+1)=yn(j) + (h/2)*(k1+k2);
end;

hold on;
plot(xn,yn,'ro')
title('exact y(t) ---   approximate y(t) ooo')
return;