MATLAB: Solution to differential equations system looks strange. What is wrong with the code

differential quationsdsolveodeprobabilities

I am writing a code to solve the differential equations. I have written the following code but my solution looks wierd. First, I get long digits numbers, which I do not get if I solve the system by hand. secondly, I get constants like C1 and C2 etc… Though I have defined the values of all parameters of functions. Please tell me what is wrong in my code:
%defined parameters 
hf12 = exp(-21.6079+0.0099).*exp(-0.0134.*t) 
hf13 = 0
hf14 = 0
hf15 = exp(-9.65573+0.01844+0.0821).*exp(0.02246*t)
hf11 = -(hf12+hf15)
hf21 = 0
hf23 = exp(-1.666-0.1116).*exp(-0.0025.*t)
hf24 = exp(-8.96236+0.07691).*exp(0.00978.*t)
hf25 = exp(-9.65573+0.08218).*exp(0.02246.*t)
hf22 = -(hf23+hf24+hf25)
   hf31 = 0
   hf32 = exp(-21.6079+0.0676+0.0099).*exp(-0.0134.*t) 
   hf34 = 0
   hf35 = exp(-9.65573-(0.11853)+0.08218).*exp(0.02246.*t)
   hf33 = -(hf32+hf35)
  hf41 = 0
hf42 = exp(-21.609-0.4176+0.0099).*exp(-0.0134.*t)
hf43 = 0
hf45 = exp(-9.65573-0.00415+0.08218).*exp(0.02246.*t)
hf44 = -(hf42+hf45)
%% solution 
syms p11(t) p22(t) p33(t) p44(t) p55(t) p12(t) p13(t) p14(t) p15(t) p21(t) p23(t) p24(t) p25(t) ...
     p31(t) p32(t) p34(t) p35(t) p41(t) p42(t) p43(t) p45(t) p51(t) p52(t) p53(t) p54(t) p55(t)
ode = [diff(p11,t) == exp(hf12+hf15),...
        diff(p22,t) == exp(hf23+hf24+hf25),...
        diff(p33,t) == exp(+hf32+hf35),...
        diff(p44,t) == exp(+hf42+hf45),...
        diff(p12,t) == p11.*hf12+p12*hf22,...
        diff(p13,t) == 0,...
        diff(p14,t) == 0,...
        diff(p15,t) == p11.*hf15+p12*hf25,...
        diff(p21,t) == 0,...
        diff(p23,t) == p22.*hf23+p23*hf33,...
        diff(p24,t) == p22.*hf24+p24*hf44,...
        diff(p25,t) == p22.*hf25+p23*hf35+p24.*hf45,...
        diff(p31,t) == 0,...
        diff(p32,t) == p32.*hf22+p33*hf32,...
        diff(p34,t) == 0,...
        diff(p35,t) == p32.*hf25+p33*hf35,...
        diff(p41,t) == 0,...
        diff(p42,t) == p42.*hf22+p44*hf42,...
        diff(p43,t) == 0,...
        diff(p45,t) == p42.*hf25+p44*hf45,...
        diff(p51,t) == 0,...
        diff(p52,t) == 0,...
        diff(p53,t) == 0,...
        diff(p54,t) == 0,...
        diff(p55,t) == 1]
    cond1 = p11(t)+p12(t)+p13(t)+p14(t)+p15(t) == 1;
    cond2 = p21(t)+p22(t)+p23(t)+p24(t)+p25(t) == 1;
    cond3 = p31(t)+p32(t)+p33(t)+p34(t)+p35(t) == 1;
    cond4 = p41(t)+p42(t)+p43(t)+p44(t)+p45(t) == 1;
    cond5 = p51(t)+p52(t)+p53(t)+p54(t)+p55(t) == 1;
S = dsolve(ode)
These equations are solutions to probabilities therefore, I am expecting the answers between 0 and 1.

Best Answer

Q = @(v) sym(v)
Q = function_handle with value:
    @(v)sym(v)

syms t real
%defined parameters 
hf12 = exp(Q(-21.6079)+Q(0.0099)).*exp(Q(-0.0134).*t) 
hf12 = 
hf13 = 0
hf13 = 0
hf14 = 0
hf14 = 0
hf15 = exp(Q(-9.65573)+Q(0.01844)+Q(0.0821)).*exp(Q(0.02246)*t)
hf15 = 
hf11 = -(hf12+hf15)
hf11 = 
hf21 = 0
hf21 = 0
hf23 = exp(Q(-1.666)-Q(0.1116)).*exp(Q(-0.0025).*t)
hf23 = 
hf24 = exp(Q(-8.96236)+Q(0.07691)).*exp(Q(0.00978).*t)
hf24 = 
hf25 = exp(Q(-9.65573)+Q(0.08218)).*exp(Q(0.02246).*t)
hf25 = 
hf22 = -(hf23+hf24+hf25)
hf22 = 
   hf31 = 0
hf31 = 0
   hf32 = exp(Q(-21.6079)+Q(0.0676)+Q(0.0099)).*exp(Q(-0.0134).*t) 
hf32 = 
   hf34 = 0
hf34 = 0
   hf35 = exp(Q(-9.65573)-Q(0.11853)+Q(0.08218)).*exp(Q(0.02246).*t)
hf35 = 
   hf33 = -(hf32+hf35)
hf33 = 
  hf41 = 0
hf41 = 0
hf42 = exp(Q(-21.609)-Q(0.4176)+Q(0.0099)).*exp(Q(-0.0134).*t)
hf42 = 
hf43 = 0
hf43 = 0
hf45 = exp(Q(-9.65573)-Q(0.00415)+Q(0.08218)).*exp(Q(0.02246).*t)
hf45 = 
hf44 = -(hf42+hf45)
hf44 = 
%% solution 
syms p11(t) p22(t) p33(t) p44(t) p55(t) p12(t) p13(t) p14(t) p15(t) p21(t) p23(t) p24(t) p25(t) ...
     p31(t) p32(t) p34(t) p35(t) p41(t) p42(t) p43(t) p45(t) p51(t) p52(t) p53(t) p54(t) p55(t)
ode = [diff(p11,t) == exp(hf12+hf15),...
        diff(p22,t) == exp(hf23+hf24+hf25),...
        diff(p33,t) == exp(+hf32+hf35),...
        diff(p44,t) == exp(+hf42+hf45),...
        diff(p12,t) == p11.*hf12+p12*hf22,...
        diff(p13,t) == 0,...
        diff(p14,t) == 0,...
        diff(p15,t) == p11.*hf15+p12*hf25,...
        diff(p21,t) == 0,...
        diff(p23,t) == p22.*hf23+p23*hf33,...
        diff(p24,t) == p22.*hf24+p24*hf44,...
        diff(p25,t) == p22.*hf25+p23*hf35+p24.*hf45,...
        diff(p31,t) == 0,...
        diff(p32,t) == p32.*hf22+p33*hf32,...
        diff(p34,t) == 0,...
        diff(p35,t) == p32.*hf25+p33*hf35,...
        diff(p41,t) == 0,...
        diff(p42,t) == p42.*hf22+p44*hf42,...
        diff(p43,t) == 0,...
        diff(p45,t) == p42.*hf25+p44*hf45,...
        diff(p51,t) == 0,...
        diff(p52,t) == 0,...
        diff(p53,t) == 0,...
        diff(p54,t) == 0,...
        diff(p55,t) == 1].';
ode = ode(t);
    cond1 = p11(t)+p12(t)+p13(t)+p14(t)+p15(t) == 1;
    cond2 = p21(t)+p22(t)+p23(t)+p24(t)+p25(t) == 1;
    cond3 = p31(t)+p32(t)+p33(t)+p34(t)+p35(t) == 1;
    cond4 = p41(t)+p42(t)+p43(t)+p44(t)+p45(t) == 1;
    cond5 = p51(t)+p52(t)+p53(t)+p54(t)+p55(t) == 1;
S = dsolve(ode)
S = struct with fields:
    p11: [1×1 sym]
    p12: [1×1 sym]
    p13: [1×1 sym]
    p14: [1×1 sym]
    p15: [1×1 sym]
    p21: [1×1 sym]
    p22: [1×1 sym]
    p23: [1×1 sym]
    p24: [1×1 sym]
    p25: [1×1 sym]
    p31: [1×1 sym]
    p32: [1×1 sym]
    p33: [1×1 sym]
    p34: [1×1 sym]
    p35: [1×1 sym]
    p41: [1×1 sym]
    p42: [1×1 sym]
    p43: [1×1 sym]
    p44: [1×1 sym]
    p45: [1×1 sym]
    p51: [1×1 sym]
    p52: [1×1 sym]
    p53: [1×1 sym]
    p54: [1×1 sym]
    p55: [1×1 sym]

S.p11
ans = 
vpa(S.p11, 5)
ans = 
Sc = struct2cell(S);
Sv = simplify(vertcat(Sc{:}), 'steps', 10)
Sv = 
symvar(Sv)
ans = 
conds = [cond1; cond2; cond3; cond4; cond5]
conds = 
temp = subs(conds, S);
%do not simplify() the equations directly, as it will test to see
%whether the left side equals the right side and will say "symfalse"
conds_subs = simplify(lhs(temp)-rhs(temp), 'steps', 10)
conds_subs = 
symvar(conds_subs)
ans = 
Related Question