Why am I unable to solve multivariate ODE boundary value problems using DSOLVE with MuPAD in Symbolic Math Toolbox 5.1 (R2008b)?

I am trying to solve the following system of second order ODEs:

 D2T1/5 - 2*DT1 = 0

 D2T2/5 - 2*DT2 + 5000 = 0

 D2T3/5 - 2*DT3 = 0

 with initial/boundary conditions:

 T1(0) = 20

 T3(1) = 100

 T1(.2) = T2(.2)

 DT1(.2) = DT2(.2)

 T2(.4) = T3(.4)

 DT2(.4) = DT3(.4)

While the following command used to work in older versions of the toolbox (with Maple), I do not receive a solution with MuPAD. When I execute

Soln = dsolve('D2T1/5 - 2*DT1 = 0', 'D2T2/5 - 2*DT2 + 5000 = 0',...

'D2T3/5 - 2*DT3 = 0',...

'T1(0) = 20', 'T3(1) = 100', 'T1(.2) = T2(.2)',...

'DT1(.2) = DT2(.2)', 'T2(.4) = T3(.4)', 'DT2(.4) = DT3(.4)');

I receive the error

 ??? Error using ==> dsolve at 145

 There are more ODEs than variables.

Upon setting a breakpoint in "dsolve.m" as advised in the related solution, I see that the error thrown by MuPAD is,

 Error: dependent variable not allowed in rhs of init conds [ode::solvesysinit]

Answer 1

[MathWorks Support Team]

This enhancement has been incorporated in Release 2009b (R2009b). For previous product releases, read below for any possible workarounds:

The ability to solve a system of ODEs with non-constant boundary conditions is not available in Symbolic Math Toolbox 5.1 (R2008b). The documentation for DSOLVE states that the initial conditions may be specified by equations like y(a)=b or Dy(a)=b where a and b are constants. Maple may be able to accept a more general type of boundary conditions, but since MuPAD is not explicitly stated to handle such specifications, it cannot be guaranteed to work with such specifications.

You may be able to work around this by first solving the system without the flux conditions to receive a solution involving constants. You can then use the SOLVE function with these conditions to solve the algebraic system of equations to determine the values of the constants and then substitute these back in the original solution.

syms t

soln = dsolve('D2T1/5 - 2*DT1 = 0', 'D2T2/5 - 2*DT2 + 5000 = 0',...

'D2T3/5 - 2*DT3 = 0', 'T1(0) = 20', 'T3(1) = 100');

eqs{1} = subs(soln.T1, t, .2) - subs(soln.T2, t, .2); % = 0

eqs{2} = subs(soln.T2, t, .4) - subs(soln.T3, t, .4);

eqs{3} = subs(diff(soln.T1,t), t, .2) - subs(diff(soln.T2,t), t, .2);

eqs{4} = subs(diff(soln.T2,t), t, .4) - subs(diff(soln.T3,t), t, .4);

constvals = solve(eqs{:});

constnames = fieldnames(constvals);

Soln.T1 = subs(soln.T1, constnames.', struct2cell(constvals).');

Soln.T2 = subs(soln.T2, constnames.', struct2cell(constvals).');

Soln.T3 = subs(soln.T3, constnames.', struct2cell(constvals).');