%type function zprime = dpdt(t,z)
function zprime = dpdt(t,z)

% for the equation dP/dt = r*(M-P)*P  we first tell them what r and M are,
% and then say what zprime is:
% r = .1;
% M = 3;
% zprime = r*(M-z)*z;

% for x'' = - x - beta x^3 + gamma cos(w t)



% for the following values of beta, gamma, and w, have three orbits
% if k=0:
%   a = -.6734177811 .1041299881 .5692877930  and all b = 0.
%
% if k = .005 then 
% a,b = (.1030367663,.01045779210)
% a,b = (.4788020648,.3291538008)
% a,b = (-.5059406275,.4291009820)
% 
% The above information is here so that I know what initial data to 
% take to make interesting phase planes.

beta = -.167;
w = .975;
g = .005;
k = 0;

x = z(1);
y = z(2);

zprime(1) = y;
zprime(2) = -x-beta*x^3+g*cos(w*t) - k*y;

% to study the van der pol plane, see equations 7.16 and 7.17
a = z(1);
b = z(2);
zprime(1) = -b/(2*w)*(w^2-1-3/4*beta*(a^2+b^2));
zprime(2) = a/(2*w)*(w^2-1-3/4*beta*(a^2+b^2)) + g/(2*w);

% for nonzero friction:
k = .005;
% the fixed points are
%  {b = .01045779210, a = .1030367663},
%  {b = .3291538009, a = .4788020649},
%  {a = -.5059406276, b = .4291009821}

zprime(1) = -4*b*w*k^2-8*b*w^3-4*k*a*w^2-4*a*k-3*k*beta*a^3-3*k*beta*a*b^2+4*k*g+8*w*b+6*beta*a^2*b*w+6*beta*b^3*w;
zprime(1) = zprime(1)/(4*k^2+16*w^2);
zprime(2) = -4*k*b*w^2+4*k^2*a*w-4*k*b-3*k*beta*a^2*b-3*k*beta*b^3+8*a*w^3-8*a*w-6*beta*a^3*w-6*beta*a*b^2*w+8*g*w;
zprime(2) = zprime(2)/(4*k^2+16*w^2);

zprime = zprime';