Invalid indexing or function definition.
Show older comments
clc;
clear all;
close all;
syms y1(t) y2(t) y3(t) y4(t) y5(t) y6(t) y7(t) y8(t)
C = 3*10^10;
KT = 0.086173324; % MeV
M0 = 1.0078; % u
M1 = 1.0078; % u
M2 = 2.014102; % u
M3 = 1.0078; % u
M4 = 2.014102; % u
M5 = 2.014102; % u
M6 = 3.016030; % u
M7 = 2.014102; % u
M8 = 3.016030; % u
M9 = 3.016030; % u
M10 = 3.016030; % u
M11 = 4.002603; % u
Z0 = 1;
Z1 = 1;
Z2 = 1;
Z3 = 1;
Z4 = 1;
Z5 = 1;
Z6 = 2;
Z7 = 1;
Z8 = 2;
Z9 = 2;
Z10 = 2;
Z11 = 2;
Eat = 0.98951013;
m01 = (M0*M1)/(M0+M1);
m23 = (M2*M3)/(M2+M3);
m45 = (M4*M5)/(M4+M5);
m67 = (M6*M7)/(M6+M7);
m89 = (M8*M9)/(M8+M9);
m1011 = (M10*M11)/(M10+M11);
cont = (931.494)./(3*10^10)^2;
cont1 = 10^-24;
sy0 = 0.1;
Q2 = (4.08*10^-15)*((y8./KT)^(-2/3))*exp(-3.381*((y8./KT)^(-1/3)))*(1+(3.82*y8./KT)+1.51*((y8./KT)^2)+0.144*(y8./KT)^3-1.14*(10^-2)*(y8./KT)^4)./(6.023*10^23);
Q4 =(2.58*10^3)*((y8./KT)^(-2/3))*exp(-3.721*((y8./KT)^(-1/3)))*(1+(3.96*y8./KT)+0.116*(y8./KT)^2)./(6.023*10^23);
Q5 = ((5.59*10^10)*(((y8./KT))^(-2/3))*exp(-12.277*((y8./KT)^(-1/3)))*(1-0.135*(y8./KT)+2.54*0.01*((y8./KT)^2)-1.29*0.001*((y8./KT)^3)))./(6.023*10^23);
D = (2.07*10^14)*((y8./KT)^(3/2))*(1-0.860*((y8./KT)^(1/2))+0.429*(y8./KT))*exp(-25.82/(y8./KT));
nB = (0.6*10^-9)*(((2*1.202)./(pi.^2))*(((10^6)*sy0.*8065.73).^3));
nP = 0.75*nB;
nN = 0.25*nB;
nEn = ((((3*1.202)./(2*pi.^2))*(((10^6)*sy0.*8065.73).^3)))+ nP;
nEp = ((((3*1.202)./(2*pi.^2))*(((10^6)*sy0.*8065.73).^3)))*(exp((-0.511)./(sy0)));
Wtn2 = 1000;
Wconst2 = ((1.239*10^-10)^3);
Wb2 = ((7.56*10^-16)*((1.239*10^-6)^3))./((1.6*10^-19)*((8.617*10^-5)^4));
WV2 = (3*10^10)*sqrt(1-(1./((y8/0.511)+1))^2);
WF2 = (((1.166*10^-23)./(8065.73*10^6))^2)*((y8)^2)*(exp(-(0.783)./y8));
WF1 = (((1.166*10^-23)./(8065.73*10^6))^2)*((y8)^2);
L1 = (y3./Wtn2)+(y5*D)-(y1*y4*WF2*WV2)-(y1*y2*WF1*WV2);
L2 = (((y4)^2)*Q2*0.5)-(y2*y3*WF1*WV2)-(y1*y2*WF1*WV2);
L3 = (y1*y4*WF2*WV2)+(y5*D)-(y2*y3*WF1*WV2)-(y3./Wtn2);
L4 = (y3./Wtn2)+(y2*y3*WF1*WV2)+(y5*D)+(((y6)^2)*Q5)-(y1*y4*WF2*WV2)-(((y4)^2)*Q2*0.5)-(y4*y5*Q4);
L5 = (((y4)^2)*Q2*0.5)-(y5*D)-(y4*y5*Q4);
L6 = (y4*y5*Q4)-(((y6)^2)*Q5*0.5);
L7 = (((y6)^2)*Q5*0.5);
S1 = (((((y4)^2)*Q2*0.5)*1.44)-((y5*D)*2.22)+((y4*y5*Q4)*5.49)+((((y6)^2)*Q5*0.5)*12.86));
S2 = (((y2*y3*WF1*WV2))+((y1*y4*WF2*WV2))+(2*(y1*y2*WF1*WV2))+((y3./Wtn2))+((((y4)^2)*Q2*0.5)))*y8;
L8 = ((Wconst2).*(S1-S2))./(4*Wb2.*((y8)^3));
ode1 = diff(y1)==L1;
ode2 = diff(y2)==L2;
ode3 = diff(y3)==L3;
ode4 = diff(y4)==L4;
ode5 = diff(y5)==L5;
ode6 = diff(y6)==L6;
ode7 = diff(y7)==L7;
ode8 = diff(y8)==L8;
odes = [ode1;ode2;ode3;ode4;ode5;ode6;ode7;ode8];
% S = dsolve(odes)
cond1 = y1(0)==nEn;
cond2 = y2(0)==nEp;
cond3 = y3(0)==nN;
cond4 = y4(0)==nP;
cond5 = y5(0)==0;
cond6 = y6(0)==0;
cond7 = y7(0)==0;
cond8 = y8(0)==0.1;
conds = [cond1;cond2;cond3;cond4;cond5;cond6;cond7;cond8];
[y1Sol(t),y2Sol(t),y3Sol(t),y4Sol(t),y5Sol(t),y6Sol(t),y7Sol(t),y8Sol(t)] = dsolve(odes,conds);
t = linspace(10^-25,10^-22,10);
loglog(t,y1Sol(t),'k-','linewidth',2)
Answers (1)
Cris LaPierre
on 4 Nov 2021
Remove the (t) from the output variables for dsolve. It is treating t as an index.
[y1Sol,y2Sol,y3Sol,y4Sol,y5Sol,y6Sol,y7Sol,y8Sol] = dsolve(odes,conds);
2 Comments
mohamed elkasemy
on 4 Nov 2021
Cris LaPierre
on 4 Nov 2021
The reason is because of the warning message that is still there (Unable to find symbolic solution). t is a 1x10 vector, but y1Sol is an empty vector.
Categories
Find more on Mathematics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!