CA_Forced_QuasiPer.m

function dxy = diffxy(t, xy)
%
%split xy into variables in our equations
%
%parameter beta = alpha/(sqrt(s16^2+c16^2))
A = xy(1);
P = xy(2);
%
% define the derivatives of these variables from equations
%
%time t is in kyr
om1=0.153249478547167;
s1=-11.2287376815124;
c1= 3.51682075211241; 
om2=0.158148666238883;
s2=-3.8249937146754;
c2=-0.761851750263805;
om3=0.117190147169570;
s3= 2.28814805956066;
c3= 1.80233702684623;
om4=0.155061775112933;
s4=-1.29770081956440;
c4=-0.635152963728496;
om5=0.217333905941751;
s5= 0.380973541305497;
c5=-1.46301711999210;
om6=0.150162587421217;
s6=1.54904176353302;
c6=-0.0883941912769817;
om7=0.211709630908568;
s7=-0.810768209286259;
c7=-0.577980646565494;
om8=0.156336369673117;
s8=-0.918358442095885;
c8=0.196083726889428;
om9=0.148350290855451;
s9=0.256895610735773;
c9=-0.524697312305024; 
om10=0.206924898030688;
s10=-0.335783913402678;
c10=-0.0194792150128644;
om11=0.212525165090383;
s11=0.267659228540196;
c11=0.128915417116900; 
om12=0.229992875969202;
s12= 0.0696189733188958;
c12= 0.0746231714061285;
om13=0.306498957094334;
s13= 0.0247349748169616;
c13= 0.0140464395340974;
om14=0.311398144786051;
s14= 0.0138353727621181;
c14= 0.0304736668840422;
om15=0.004899187691716;
s15=-0.160479848721994;
c15=0.0594077968934257;
om16=0.264933601588513;
s16=-15.5490493322904;     
c16=-9.70406287110532;
om17=0.280151350350945;
s17=15.4319556361701;
c17= 4.75247271131525;
om18=0.331110950251899;
s18= 9.0992249352734;
c18=-10.6115244887390;  
om19=0.328024059125949;
s19=-7.87065384013669;
c19= 6.61544246063503;
om20=0.326211762560183;
s20= 0.813786144754451;
c20=-4.52641408099246;
om21=0.269742342439881;
s21= 0.0690448504314857;
c21=-3.31639260969558;
om22=0.332923246817665;
s22=1.44050770785967;
c22= 1.06339286050120;
om23=0.371638925683567;
s23= 0.925324276580528;
c23=-1.02066758672154;
om24=0.275366617473065;
s24= 0.997628846513796;
c24=-0.362906496840039;  
om25=0.323124871434233;
s25=-0.378637986107629;
c25=0.527217891742183;
om26=0.259396912994958;
s26= 0.339477750517033;
c26=-0.560509461538342;
om27=0.324937167999999;
s27=-0.576082669762308; 
c27=1.18669572739338;
om28=0.334197841377850;
s28= 0.346906064369828;
c28=-0.648189701487285;
om29=0.274551083291250;
s29=-0.441772417569753;
c29= 0.289576210423804;
om30=0.418183080135680;
s30=-0.0184884064645011; 
c30=0.109632390175297;
om31=0.111684123041346;
s31=-0.428006728186239; 
c31=0.357006342316690;
om32=0.433400828898112;
s32=-0.0049199219454561;
c32=-0.106148873639336;
om33=0.126901871803777;
s33= 0.257509918217341;
c33=-0.377639794223366; 
om34=0.336010137943616;
s34=-0.421809264016129;
c34=0.324327509437558;
om35=0.177861471704732;
s35=-0.161827722328271;
c35=-0.362683869407858;
F=s1*sin(om1*t)+c1*cos(om1*t)+s2*sin(om2*t)+c2*cos(om2*t)+s3*sin(om3*t)+c3*cos(om3*t)+s4*sin(om4*t)+c4*cos(om4*t)+s5*sin(om5*t)+c5*cos(om5*t)+s6*sin(om6*t)+c6*cos(om6*t)+s7*sin(om7*t)+c7*cos(om7*t)+s8*sin(om8*t)+c8*cos(om8*t)+s9*sin(om9*t)+c9*cos(om9*t)+s10*sin(om10*t)+c10*cos(om10*t)+s11*sin(om11*t)+c11*cos(om11*t)+s12*sin(om12*t)+c12*cos(om12*t)+s13*sin(om13*t)+c13*cos(om13*t)+s14*sin(om14*t)+c14*cos(om14*t)+s15*sin(om15*t)+c15*cos(om15*t)+s16*sin(om16*t)+c16*cos(om16*t)+s17*sin(om17*t)+c17*cos(om17*t)+s18*sin(om18*t)+c18*cos(om18*t)+s19*sin(om19*t)+c19*cos(om19*t)+s20*sin(om20*t)+c20*cos(om20*t)+s21*sin(om21*t)+c21*cos(om21*t)+s22*sin(om22*t)+c22*cos(om22*t)+s23*sin(om23*t)+c23*cos(om23*t)+s24*sin(om24*t)+c24*cos(om24*t)+s25*sin(om25*t)+c25*cos(om25*t)+s26*sin(om26*t)+c26*cos(om26*t)+s27*sin(om27*t)+c27*cos(om27*t)+s28*sin(om28*t)+c28*cos(om28*t)+s29*sin(om29*t)+c29*cos(om29*t)+s30*sin(om30*t)+c30*cos(om30*t)+s31*sin(om31*t)+c31*cos(om31*t)+s32*sin(om32*t)+c32*cos(om32*t)+s33*sin(om33*t)+c33*cos(om33*t)+s34*sin(om34*t)+c34*cos(om34*t)+s35*sin(om35*t)+c35*cos(om35*t); 
k=50;
S=100;
I=0.004;
beta=0.00025;

Adot = I-k*(1+beta*F)*A*exp(P);
Pdot = k*(1+beta*F)*A-S;
%i
%return the derivatives in dxy in the right order
%
dxy = [Adot; Pdot];