I have solved the differential equation given in the attachments using ode45 now when i am trying to do it in simulink its giving me error of singularity and that too the response is not similar. Can someone identify the error that I am doing please?

The simulink file is attached in as a pdf. The code for function file and script file I am adding below.

function [ dy ] = rigid(t,y)

%data

cl=8;

cm=6;

ec=2.9442*10^3;

ed=2.9442*10^3;

em=7.4478*10^4;

el=1.2550*10^5;

ep=1.8283*10^4;

f0=.58;

m=100.12;

r=8.314;

t=335;

zc=3.8223*10^10;

zd=3.1457*10^11;

zm=1.0067*10^15;

zl=3.7920*10^18;

zp=1.7700*10^9;

v=0.1;

f=1;

global q;

%% where po here is

po=(2*f0*y(2)*zl*exp(-el/(r*t))/(zd*exp(-ed/(r*t))+zc*exp(-ec/(r*t))))^0.5;

%% here goes the code

%y(1)=Cm; y(2)=C1; y(3)=D.; y(4)=D1;

dy=zeros(4,1);

dy(1)=-(zp*exp(-ep/(r*t))+zm*exp(-em/(r*t)))*y(1)*po-f*y(1)/v+f*cm/v;

dy(2)=-zl*exp(-el/(r*t))*2.2433*10^-1-(f*y(2)/v)+q*cl/v;

dy(3)=(0.5*zc*exp(-ec/(r*t))+zd*exp(-ed/(r*t)))*po^2+zm*exp(-em/(r*t))*y(1)*po-f*y(3)/v;

dy(4)=m*(zp*exp(-ep/(r*t))+zm*exp(-em/(r*t)))*y(1)*po-f*y(4)/v;

end

%% And here is the script file for running

clc;

clear all;

global q;

%% data available steady state

x1_ss=5.374813005334617;

x2_ss=.224330229793841;

x3_ss=.003130435224154;

x4_ss=62.593721905898185;

x_ss=[x1_ss;x2_ss;x3_ss;x4_ss];

f1min=0.003;

f1max=0.06;

%% Iteration

for i=1:1

f1(i)=f1min+rand(1,1)*(f1max-f1min);

q=f1(i);

j=1*i-1;

l=j+1;

[t,y] = ode45(@rigid,[j:.1:l],x_ss);

[k,p]=size(t);

T(1:k,i)=t;

r=ones(k,i);

g(1:k,i)=f1(i)*r(:,i);

NAMW(1:k,i)=y(1:k,4)./(y(1:k,3));

x_ss=[y(k,1);y(k,2);y(k,3);y(k,4)];

end

out1=NAMW(:);

out2=T(:);

out3=g(:);

%% plots

figure(1)

plot(out2,out1);

title('training sample');

ylabel('NAMW');

xlabel('time');

figure(2)

plot(out2,out3);

title('training sample');

ylabel('F1');

xlabel('time');