function test2
% Solving damping and resonance
clear all
%initial conditions
 t0= 0; x0 = 0;  v0 = 0; 

%%parameters
m = 1;  %value of mass
k = 1;  %value of spring constant
g = 10; %gravity
c = .05;  %drag constant
h = .001;   % step size = delta t
omega= 1.02;  % forcing frequency
L0   = 0;   % rest length of spring
d     = 2;  % forcing amplitude
t0 = 0;  t1 = 50*pi;  % interval to solve odes


n= (t1-t0)/h;  %number of time steps
n = floor(n);
% initialize
tnew=t0; xnew = x0;  vnew=v0; 

tarray= zeros(1,n+1);xarray=zeros(1,n+1); varray=zeros(1,n+1);
tarray(1)=t0;      xarray(1) = x0;    varray(1) = v0;

%%%%%%%%
for i = 2:n+1
x = xnew; v = vnew; t = tnew;
F=  rhs(x,v,t,k,c,d,L0,m,g, omega);

a = F/m;          % the laws of mechanics
xnew = x + v*h;   % kinematics
vnew = v + a*h;   % kinematics
tnew = t + h;

tarray(i) = tnew;
xarray(i) = xnew;
varray(i) = vnew;
end

figure(1)
plot(tarray, xarray, 'r')
xlabel('time'); ylabel('position');
title('Andy''s plot','fontsize', 20, 'color','b')

%figure(2)
%plot(xarray,varray, 'b')
%axis('square'); axis('equal')
end


function  F=rhs(x,v,t,k,c,d,L0,m,g, omega)
F = -k*(x-d*sin(omega*t) - L0) - c*v + m*g;
end