function ode_solver()
% model driver for ode 45
% define system parameters
 m = 1;
 c = 1;
 k = 1;
 w = 2.5;
%define range of integration
 t0 = 0; tfinal = 32*pi/w;
 tspan = [t0 tfinal];
% hold on % use this to graph several solutions on same plot
 options = []; % will need to alter to increase precision etc.

 x0 = 1; y0 = 0;
 [t z] = ode45(@pendulum, tspan, [x0 y0],options, m, c, k, w)
 N = length(t); %querry length of t-vector
 x = z(:,1); y = z(:,2); % grab values for x, y; use more variables if
                         % higher dimensional system
 x(N) % to see end-value of x
 plot(x,y) % plot y vs. x; change to plot(t,x) etc. if other plots desired
 axis equal;

 function zp = pendulum(t, z, m, c, k, w)
 x = z(1); y = z(2); % for cosmetic reasons: read z(1), z(2) into x, y etc
% treat t, z; x, y, ...  as vector variables: Use ./  .^  .*  etc.
% 1st equation in system:
 xp = y;
% 2nd equation in system:
 yp = -(c/m)*y - (k/m)*sin(x) + (1/m)*cos(w*t);
% construct the zp vector from its components
 zp = [xp; yp];
