I am attempting to solve a system of second order ODEs. I place conditions on the solutions and use the solve command to figure the correct constants for the general solutions of the ODEs; however, the conditions do not appear to hold after I substitute the constants back into the general solutions. Any help would be greatly appreciated. Here's the code and an explanation:

First some constants

> A := 1; B := 9/10;

> j := 1-1/B;

This is our homogeneous odes. I will give the general solutions of the inhomogeneous system momentarily

> eqnv1 := diff(v1(x), `$`(x, 2)) = (1-1/(j+1))*v1(x)+v2(x)/(j+1);

> eqnv2 := diff(v2(x), `$`(x, 2)) = -v1(x)/(A*(j+1))+(B/A+1/(A*(j+1)))*v2(x);

Next we get the general solution of this sytem of odes.

> soln := dsolve([eqnv1, eqnv2])

Next we have our solutions of the inhomogeneous problem1. Basically solution v1neg, v2neg on [0,xi] and v1pos, v2pos on [xi,1]. We will assume v1,v2 are C^1 across xi; however, the location of xi is not known at this time so they must remain split.

> v1neg := op([1, 2], soln)-1;

> v2neg := op([2, 2], soln)-1/B;

> v1pos := op([1, 2], soln)+1;

> v2pos := op([2, 2], soln)+1/B;

There's probably a better way to do this, but I relabeled the constants:

> v1negc := subs([_C1 = a[1], _C2 = a[2], _C3 = a[3], _C4 = a[4]], v1neg);

> v2negc := subs([_C1 = a[1], _C2 = a[2], _C3 = a[3], _C4 = a[4]], v2neg);

>

> v1posc := subs([_C1 = a[5], _C2 = a[6], _C3 = a[7], _C4 = a[8]], v1pos);

> v2posc := subs([_C1 = a[5], _C2 = a[6], _C3 = a[7], _C4 = a[8]], v2pos);

Next we have eight conditions the solutions must satisfy. Namely v1, v2 are C^1 across xi and v1',v2' are 0 at {0,1}.

> syscon1 := subs(x = xi, v1negc) = subs(x = xi, v1posc);

> syscon2 := subs(x = xi, v2negc) = subs(x = xi, v2posc);

> syscon3 := subs(x = xi, diff(v1negc, x)) = subs(x = xi, diff(v1posc, x));

> syscon4 := subs(x = xi, diff(v2negc, x)) = subs(x = xi, diff(v2posc, x));

> syscon5 := subs(x = 0, diff(v1negc, x)) = 0;

> syscon6 := subs(x = 0, diff(v2negc, x)) = 0;

> syscon7 := subs(x = 1, diff(v1posc, x)) = 0;

> syscon8 := subs(x = 1, diff(v2posc, x)) = 0;

We solve to get the constants for the solutions.

> constants := simplify(evalf(solve({syscon1, syscon2, syscon3, syscon4, syscon5, syscon6, syscon7, syscon8}, {a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]})));

>

We substitute the values for the constants.

> a[1] := op([1, 2], constants); a[2] := op([2, 2], constants); a[3] := op([3, 2], constants); a[4] := op([4, 2], constants); a[5] := op([5, 2], constants); a[6] := op([6, 2], constants); a[7] := op([7, 2], constants); a[8] := op([8, 2], constants);

Lastly we try to verify that the conditions from earlier hold:

> evalf(subs(xi = .2, subs(x = xi, v1negc-v1posc)));

-1.7597825261536669519

> evalf(subs(xi = .2, subs(x = xi, v2negc-v2posc)));

-1.8936659961101033997

> evalf(subs([x = 0, xi = .2], diff(v1negc, x)));

-0.38633519704430619686

They should hold for any xi, but they don't appear to. All of these should be 0. For a large xi, the numbers get very large so I was thinking perhaps roundoff error, but even when I do an exact solution and then evalf just at the end, I still have large error so I'm not sure what the problem is. Sorry for the long question. Thanks so much for the help.