solving system of ode using matlab -

i have 9 equations time dependent coefficient g

% m file function dy =tarak(t,y) g= 3.16; g =  0.1*exp(-((t-200)/90).^2); dy=zeros(9,1); dy(1)=-2*2*y(1)+2*g*y(5)+2*g*y(7); dy(2)=2*y(1)-2*g*y(5); dy(3)=2*y(1)-2*g*y(7); dy(4)=-2*y(4)+g*y(9); dy(5)=-2*y(5)+g*(y(2)-y(1))+g*y(8); dy(6)=-2*y(6)-g*y(9); dy(7)=-2*y(7)+g*(y(3)-y(1))+g*y(8); dy(8)=-g*y(7)-g*y(5); dy(9)=g*y(6)-g*y(4); 

then in command window:

[t,y] = ode45(@tarak,[0 ,500],[0 0 1 0 0 0 0 0 0]) 

where coefficient g = 3.16 , g = 0.1*exp(-((t-200)/90).^2) time dependent coefficient , time t = 0:500; initial condition [0 0 1 0 0 0 0 0 0].

i'm getting wrong negative values output y(1), y(2). can pls try solve above eqns ode45 can compare results.

with simple application of rk4 picture

nicely positive, 1 strange initial jump in y(1) component. note scale, on whole y(1) rather small. seems system stiff @ point, rk45 might have problems, implicit runge-kutta method better.

and zoom of initial oscillations

---

python code

import numpy np import matplotlib.pyplot plt math import exp  def dydt(t,y):     dy = np.array(y);      g = 3.16;     g = 0.1*exp(-((t-200)/90)**2);      dy[0]=-2*2*y[0]+2*g*y[4]+2*g*y[6];     dy[1]=   2*y[0]-2*g*y[4];     dy[2]=   2*y[0]-2*g*y[6];     dy[3]=  -2*y[3]+  g*y[8];     dy[4]=  -2*y[4]+  g*(y[1]-y[0])+g*y[7];     dy[5]=  -2*y[5]-  g*y[8];     dy[6]=  -2*y[6]+  g*(y[2]-y[0])+g*y[7];     dy[7]=  -g*y[6]-  g*y[4];     dy[8]=   g*y[5]-  g*y[3];     return dy;  def rk4step(f,x,y,h):     k1=f(x      , y         )     k2=f(x+0.5*h, y+0.5*h*k1)     k3=f(x+0.5*h, y+0.5*h*k2)     k4=f(x+    h, y+    h*k3)     return (k1+2*(k2+k3)+k4)/6.0   t= np.linspace(0,500,200+1); dt = t[1]-t[0]; y0=np.array([0, 0, 1, 0, 0, 0, 0, 0, 0]);  y = [y0]  t0 in t[0:-1]:     n=200;     h = dt/n;     in range(n):         y0 = y0 + h*rk4step(dydt,t0+i*h,y0,h);     y.append(y0);  y = np.array(y);  plt.subplot(121); plt.title("y(1)") plt.plot(t,y[:,0],"b.--") plt.subplot(122); plt.title("y(2)") plt.plot(t,y[:,1],"b-..") plt.show()