%Huen's method - second-order Runge Kutta
%Sove stiff ODE y'= -100*y+100*t+101, initial value (should be close to 1)
%time step and duration of interation, are input paramethers.
% plot solution y(t)

function [] = huen1 (y0,h,time)
c=y0-1; % coresponcing constant in the general solution
xp= linspace(0,time,100); %always take 100 points
yp= (1+xp+c*exp(-100*xp)); % 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= -100*yn(j)+100*xn(j)+101;
k2= -100*(yn(j)-h*k1)+100*xn(j+1)+101; 
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;