做一个课题，需要计算很多微分方程，之前一直用MATLAB，前面的方程还行，后面的方程是刚性方程所以解有问题，所以看了下科普，MATLAB算是比较低端的软件了，所以转向mathematica试试。
结果出乎意料，之前MATLAB能解的方程，mathematica居然不能解
MATLAB解出来的结果带有wrightOmega项来表示，这好像是MATLAB自创的一种表示，意思是方程是超越方程，没有解析代数表达式，但是可以代入具体数得到数值解。然而mathematica直接报“病态矩阵”，啥都解不出来，也不能带具体数值进行计算...说好的MATLAB才是最低端的入门级软件呢，MATLAB都能解出来的mathematica居然搞不定，估计后面的MATLAB解不了的就更不行了。总之再试试maple

贴具体的代码！

matlab不低端，mma也需要添加一些选项才能得到解。

c = 299792458*10^2(*光速，单位cm/s*)
G = 6.67259*10^-8(*gravitational constant，引力常数，单位cm^3/g*s^2*)
Msun = 1.9891*10^33(*Subscript[M, \[CircleDot]]，太阳质量，单位g*)
Mbi = 1.4*Msun(*Subscript[M, b,i]，(29)式下面*)
\[Eta] = 0.01(*\[Eta]，fig1的不同情况*)
t0 = 3000(*Subscript[t, 0]，fig1的不同情况*)
R = 11.5*10^5(*单位cm，R=11.5km，(29)式下面*)
Medot[t_] :=
10^-3*\[Eta]*t^(1/2)*Msun(*(13)式，Subscript[Overscript[M, .], early]*)
Mldot[t_] :=
10^-3*\[Eta]*t0^(13/6)*t^(-5/3)*
Msun(*(14)式，Subscript[Overscript[M, .], late]*)
Mdot[t_] := (1/Medot[t] + 1/Mldot[t])^-1(*(12)式，Overscript[M, .]*)
M[t_] := Mb[t]*(1 + 3*G*Mb[t]/(5*R*c^2))^-1(*(16)式，M*)
Mbdot[t_] :=
Mdot[t]/((1 + 3/5*G*Mb[t]/(R*c^2))^-1 -
Mb[t]*3/5*G/(R*c^2)*(1 + 3/5*G*Mb[t]/(R*c^2))^-2)
DSolve[{Mb'[t] == Mbdot[t], Mb[0] == Mbi}, Mb[t], t]
既然要我贴，那就贴吧。不知道直接复制粘贴的格式会不会丢失信息。倒是希望是我自己有不知道的额外选项能解这个方程

如果Maple也让你失望了，那你怎么办？你确定你会用Mathematica吗，刚入门就想指点江山了。

Medot[t_]:=(\[Eta] Sqrt[t] Msun)/10^3;
Mldot[t_]:=(\[Eta] t0^(13/6) t^(-5/3) Msun)/10^3;
Mdot[t_]:=1/(1/Medot[t]+1/Mldot[t]);
M[t_]:=Mb[t]/(1+(3 G Mb[t])/(5 R c^2));
Mbdot[t_]:=Mdot[t]/(1/(1+(3 G Mb[t])/(5 (R c^2)))-(Mb[t] 3 G)/(5 (R c^2) (1+(3 G Mb[t])/(5 (R c^2)))^2));DSolve[{Mb'[t]==Mbdot[t],Mb[0]==Mbi}//Simplify,Mb[t],t]

先求解再赋值解析解都能算出来。。。作为高等生物请把软件这类非生物东西熟悉了再喷，它也很无辜，不管Mma还是matlab或maple。

Clear["Global`*"]
c = 299792458*10^2(*光速，单位cm/s*)
G = 6.67259*10^-8(*gravitational constant，引力常数，单位cm^3/g*s^2*)
Msun = 1.9891*10^33(*Subscript[M, \[CircleDot]]，太阳质量，单位g*)
Itilder = 0.283(*Overscript[I, ~]，(7)式下面*)
Jtilder = 1.81*10^-2(*Overscript[J, ~]，(23)式下面*)
Mbi = 1.4*Msun(*Subscript[M, b,i]，(29)式下面*)
\[Eta] = 0.01(*\[Eta]，fig1的不同情况*)
t0 = 3000(*Subscript[t, 0]，fig1的不同情况*)
B = 2*10^14(*(29)式下面，单位G*)
R115 = 1(*Subscript[R, 11.5]*)
R = 11.5*10^5(*单位cm，R=11.5km，(29)式下面*)
M14 = 1(*Subscript[M, 1.4]*)
B15 = 1(*Subscript[B, 15]*)
Bt14 = 1(*Subscript[B, t,14]*)
T9 = 1(*Subscript[T, 9]*)
u = B*R^3(*\[Mu]，磁偶极矩，(18)式下面*)
Medot[t_] :=
10^-3*\[Eta]*t^(1/2)*Msun(*(13)式，Subscript[Overscript[M, .], early]*)
Mldot[t_] :=
10^-3*\[Eta]*t0^(13/6)*t^(-5/3)*
Msun(*(14)式，Subscript[Overscript[M, .], late]*)
Mdot[t_] := (1/Medot[t] + 1/Mldot[t])^-1(*(12)式，Overscript[M, .]*)
M[t_] := Mb[t]*(1 + (3*G*Mb[t])/(5*R*c^2))^-1(*(16)式，M*)
Mbdot[t_] :=
Mdot[t]/((1 + 3/5*G*Mb[t]/(R*c^2))^-1 -
Mb[t]*3/5*G/(R*c^2)*(1 + 3/5*G*Mb[t]/(R*c^2))^-2)
II[t_] := Itilder*M[t]*R^2(*I，(6)式下面*)
v3[t_] := \[CapitalOmega][t]/(2*Pi*10^3)(*Subscript[\[Nu], 3]，(2)式下面*)
CC[t_] := 2.3*10^-4*Bt14^2*M14^-1*v3[t]^-2(*C，改写为CC(7)式下面*)
rm[t_] := (u^4/(G*M[t]*Mbdot[t]^2))^(
1/7)(*Subscript[r, m]，(18)式，单位km*)
\[Omega][t_] := \[CapitalOmega][t]/Sqrt[G*M/rm[t]^3](*\[Omega]，(19)式*)
n\[Omega][t_] := 1 - \[Omega][t](*n(\[Omega])，(21)式下面*)
rc[t_] := 16.5*M14^(1/3)*v3[t]^(-2/3)(*Subscript[r, c]，(17)式，单位km*)
tsv = 1.4*10^8*M14*R115^-1*T9^(5/3)(*Subscript[t, sv]，(3)式，单位s*)
tgw[t_] := -24*R115^-4*M14^-1*v3[t]^-6(*Subscript[t, gw]，(2)式，单位s*)
tBt[t_] :=
5.8*R115^-1*M14*B15^-1*Bt14^-1*v3[t](*Subscript[t, B,t]，(5)式，单位s*)
tBgw[t_] :=
3.8*10^14*M14^-1*R115^-2*Bt14^-4*
v3[t]^-4(*Subscript[t, B,gw]，(8)式，单位s*)
Nacc[t_] := n\[Omega][t]*u^2/(rm[t]^3)(*(20)式*)
tdip[t_] := -1.4*10^3*B15^-2*v3[t]^-2(*(10)式，单位s*)
tacc[t_] := II[t]*\[CapitalOmega][t]/Nacc[t](*(21)式，单位s*)
tbv[t_] := 1/(1.4*10^-10*R115^5*M14^-1*T9^6*
v3[t]^2*(1 + 284.5*R115^4*v3[t]^4*T9^-2*\[Alpha][t]^2 +
3.16*10^4*R115^8*v3[t]^8*T9^-4*\[Alpha][t]^4 +
1.08*10^6*R115^12*v3[t]^12*T9^-6*\[Alpha][t]^6))(*(4)式*)
Aplus[t_] := 1 + (3*\[Alpha][t]^2*Jtilder)/(2*Itilder)(*(28)式下面*)
Aminus[t_] := 1 - (3*\[Alpha][t]^2*Jtilder)/(2*Itilder)(*(28)式下面*)
d\[Alpha][
t_] := \[Alpha][
t]*((Mbdot[t]*Aminus[t])/(2*M[t]*Aplus[t]) -
Aminus[t]/Aplus[t]*(1/tsv + 1/tbv[t] + 1/tBt[t]) -
1/Aplus[t]*(1/tdip[t] + 1/tacc[t] + 1/tBgw[t]) - 1/
tgw[t])(*\[Alpha]'，(28)式*)
d\[CapitalOmega][
t_] := \[CapitalOmega][
t]*(-(Mbdot[t]/(Aplus[t]*M[t])) +
1/Aplus[t]*(1/tdip[t] + 1/tacc[t] + 1/tBgw[t]) - (
3*\[Alpha][t]^2*Jtilder)/(
Itilder*Aplus[t])*(1/tsv + 1/tbv[t] + 1/
tBt[t]))(*\[CapitalOmega]，(27)式*)
dBt[t_] := (4/(3*Pi))^(1/2)*
B*\[Alpha][t]^2*\[CapitalOmega][t](*Subscript[B, t]，(6)式*)
NDSolve[{Mb'[t] == Mbdot[t], \[Alpha]'[t] ==
d\[Alpha][t], \[CapitalOmega]'[t] == d\[CapitalOmega][t],
Bt'[t] == dBt[t],
Mb[0] == Mbi, \[Alpha][0] == 10^-8, \[CapitalOmega][0] == (2*Pi)/(
3*10^3), Bt[0] == 100}, {Mb[t], \[Alpha][t], \[CapitalOmega][t],
Bt[t]}, {t, 0, 10000}]
还是不行啊，前面的问题是遇到无穷大，忽略了也能得到图形，但是加上后面的方程组成微分方程组就完全不行了

\[Omega][t_] := \[CapitalOmega][t]/Sqrt[G*M[t]/rm[t]^3](*\[Omega]，(19)式*)
...
expr = Simplify[{Mb'[t] == Mbdot[t], \[Alpha]'[t] ==
d\[Alpha][t], \[CapitalOmega]'[t] == d\[CapitalOmega][t],
Bt'[t] == dBt[t],
Mb[0] == Mbi, \[Alpha][0] ==
10^-8, \[CapitalOmega][0] == (2*Pi)/(3*10^3), Bt[0] == 100}];
NDSolve[expr, {Mb[t], \[Alpha][t], \[CapitalOmega][t], Bt[t]}, {t, 0.,
10000}] // Quiet
(*见不得警示不如不看*)

一看代码就知道其实是自己不会用。

建议楼主看一下
Astrophysics Through Computation: With Mathematica® Support
这本书

MATLAB 低端吗？