solution to differential equations system looks strange. What is wrong with my code?

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.

[susman]

where, t can be any time period like t = 20

Answer 1

[Walter Roberson]

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 = 

[Walter Roberson]

So the conditions involve all 25 constants. The constants exist because you define 25 differential equations of 25 functions, and each function needs a boundary condition to fully define it.

Your conditions do put constraints on the boundary conditions. With 5 equations, you could hope to solve for 5 of the constants. Which 5? Well, rememeber that they are constant -- so whatever solution is found must not involve indefinite integrals (but could potentially involve integrals over fixed boundaries, as those become constants.) Also, remember that you cannot algebraically solve for a variable that appears inside an integral (unless you can move it out of the integral.) It is not obvious to me that there is any hope in this situation.

[Walter Roberson]

Notice for example that the final condition simplifies to C1 + C2 + C3 + C4 + C5 - t = 1 . That condition is inconsistent with any of C1, C2, C3, C4, C5 actually being constant.

[susman]

Dear Walter,

Thank you for a detailed explaination. I was assuming that it would not be that complicated.

What do you suggest, how should I get the solution of it? My aim is to get the transition probabillity matrix for Markov Chain. The information I have include a Generator Matrix (transition intensity matrix) for five states,, which evolves in each period, like

G(t) = [hf11, hf12, hf13, hf14, hf15; hf21, hf22, hf23, hf24, hf25; hf31, hf32, hf33, hf34, hf35; hf41, hf42, hf43, hf44, hf45; hf51 hf52, hf53, hf54, hf55]

where each element "hf.." I have defined in my code varies with time.

I solved this equation using Kolmogorov equations. I determined the differential equations for each possible transition and then made the above code.

What do you say, its a problem of my methodology or coding?

[Walter Roberson]

Your first four p5* have derivative 0 so they must be constant. Your p55 derivative is 1, so the function must be 1*t plus a constant. The five together total 1. But the five added are a series of constants plus t. That cannot be constant unless t is restricted to 0. Therefore your equations are wrong.

[Walter Roberson]

Notice that you have exp() of terms that are already exp(). See sigma 19 in the above "where" list to see it come in as a factor. You can be fairly sure that you will not be able to find a closed form solution for integrals or ode that include such terms.