Modify Euler's Method to Heun's method

Hey all i have coded Eulers method, however i now need to modify it to include Heun's method this is what i have so far...

% % This is the work of Dr. Gretarson
% % Damped Simple Harmonic Oscillator solved by Euler's Method
% %
% % This routine implements a numerical (Euler's Method) solution to the
% % damped simple harmonic oscillator.
% % Physical Parameters of the Oscillator
% w0=10; beta=1;
% w=sqrt(w0^2-(beta/2)^2);
% phi=atan(-1/20);
% x0=1/cos(phi);
% % Parameters for the Integration
% h = 0.003; % Time step in seconds
% duration=15; % Duration of solution (seconds)
% Nstep=round(duration/h); % Resulting number of steps needed.
% t=0:h:Nstep*h; % Resulting time vector (seconds)
% % Reserve some memory for the solution vector (these lines are optional but
% % do speed up the code.)
% x=zeros(size(t)); % x is position (meters)
% y=zeros(size(t)); % y is velocity (meters)
% % Set the initial conditions to be the first elements of x and y.
% x(1)=x0; % Initial position (meters)
% y(1)=0; % Initial velocity (meters/second)
% % Iterate to find solution with Euler's method
% for n=2:Nstep
% 
%  f=y(n-1); % f(x,y,t) is dx/dt
%  x(n)=x(n-1)+f*h; % Use the last position to find the new position
%  g=-beta*y(n-1)-w0^2*x(n-1); % g(x,y,t) is dy/dt.
%  y(n)=y(n-1)+g*h; % Use the last velocity to find the new velocity
%  
%  slopery = (y(n)+h*(x(n)));
%  
%  sloperx = x(n)+h; 
%  
%  ideal = (1/2)*(sloperx + slopery);
%  
% 
% end
% % The exact solution (solved on paper) is
% xexact=x(1)*exp(-beta/2*t).*cos(w*t+phi);
% yexact=x(1)*(-beta/2*exp(-beta/2*t).*cos(w*t+phi)-w*exp(-beta/2*t).*sin(w*t+phi));
% GDE=max(xexact-x);
% % Compare the exact solution and the numerical solution
% figure(1);
% plot(t,x,'bo',t,xexact,'r','linewidth',1,'markersize',2); hold on;
% title(['GDE = ',num2str(GDE)]);
% grid on; hold off;
% % =======================================================================
% figure(2)
% plot(t,ideal,'bo',t,xexact,'r')

[Geoff Hayes]

Rob - what is your question and why is all of your code commented out?

[Rob Mullins]

So I think I got it... can someone review my code and make sure i am currect. Thank you.

% % Damped Simple Harmonic Oscillator solved by Euler's Method
% %
% % This routine implements a numerical (Euler's Method) solution to the
% % damped simple harmonic oscillator.
% % Physical Parameters of the Oscillator
% w0=10; beta=1;
% w=sqrt(w0^2-(beta/2)^2);
% phi=atan(-1/20);
% x0=1/cos(phi);
% % Parameters for the Integration
% h = 0.01; % Time step in seconds
% duration=15; % Duration of solution (seconds)
% Nstep=round(duration/h); % Resulting number of steps needed.
% t=0:h:Nstep*h; % Resulting time vector (seconds)
% % Reserve some memory for the solution vector (these lines are optional but
% % do speed up the code.)
% x=zeros(size(t)); % x is position (meters)
% y=zeros(size(t)); % y is velocity (meters)
% ideal = zeros(size(t));
% % Set the initial conditions to be the first elements of x and y.
% x(1)=x0; % Initial position (meters)
% y(1)=0; % Initial velocity (meters/second)
% % Iterate to find solution with Euler's method
% for n=2:Nstep
% 
%  f=y(n-1); % f(x,y,t) is dx/dt
%  x(n)=x(n-1)+f*h; % Use the last position to find the new position
%  g=-beta*y(n-1)-w0^2*x(n-1); % g(x,y,t) is dy/dt.
%  y(n)=y(n-1)+g*h; % Use the last velocity to find the new velocity
%  
%  slopery(n) = (y(n)+h*(x(n)));
%  
%  sloperx(n) = x(n)+h; 
%  
%  ideal(n) = (1/2)*(sloperx(n) + slopery(n));
%  
% 
% end
% % The exact solution (solved on paper) is
% xexact=x(1)*exp(-beta/2*t).*cos(w*t+phi);
% yexact=x(1)*(-beta/2*exp(-beta/2*t).*cos(w*t+phi)-w*exp(-beta/2*t).*sin(w*t+phi));
% GDE=max(xexact-x);
% % Compare the exact solution and the numerical solution
% figure(1);
% plot(t,x,'bo',t,xexact,'r','linewidth',1,'markersize',2); hold on;
% title(['GDE = ',num2str(GDE)]);
% grid on; hold off;
% % =======================================================================
% figure(2)
% plot(t,ideal,'bo',t,xexact,'r')

I have the difference in values and Eulers method is more accurate as of now... I go the difference for Euler as 1.0012 and for Heun's i got 1.0125.. I think thhis makes since as Eulers method only fails as the function is concave down and for very high digits...Am i correct?

[Rob Mullins]

Geoff, Sorry I must have not completed my thought process. My question was how do I implement the Huens method. However as seen above i might have got it but my error values are off. I dont think im getting the correct coords.