Thanks everyone for trying to help. I ended up just writing the code in sage and it works so I no longer require any assistance.

I'm interested in a discontinuous solution for modeling purposes. 

OneDiscontTwoPiecewise.mw

This is my program. The system is given by 

-v + ju + f(u) = 0

-v_1''+v_1 = u

-Av_2''+Bv_2=u

v_1'(0)=v_2'(0)=v_1'(1)=v_2'(1)

where

f(u) = u+1 for u < 0 and u-1 for u>0

v = v_1-v_2

u is assumed to have one discontinuity and v = v_1-v_2. Also v_1, v_2 are C^1. Basically you are solving the system on two pieces, one on each side of the discontinuity of u. This produces solutions for v_1, v_2 on each side of the discontinuity with four arbitrary constants (eight constants in total). Furthermore, the requirement that v_1,v_2 be C^1 or just strong solutions of the system along with the given boundary conditions gives us eight conditions for v_1,v_2 to satisfy as so we should be able to solve for the constants explicitly. Maple does find the constants; however, when substituting them in, the given solutions no longer have the correct characteristics, ie v_1'(0) is not 0.  

I do not expect u to be unique without further information. It will be determined by the conditions on the boundary and the number of discontinuities. If I know where the discontinuity is, I can uniquely solve for it on each interval of continuity. The set up of the problem isn't really what is bothering me though. It's the fact that I use solve, it finds constants where the conditions supposively hold, but then I check if the conditions hold with the solved for constants and they do not.  

I'm trying to solve the following system:

-v + ju + f(u) = 0

-v_1''+v_1 = u

-Av_2''+Bv_2=u

v_1'(0)=v_2'(0)=v_1'(1)=v_2'(1)

where

f(u) = u+1 for u < 0 and u-1 for u>0

v = v_1-v_2

u is allowed to be discontinuous and in the program I am assuming one discontinuity at xi. I would like to eventually solve the problem for arbitrarily many discontinuities in u. A,B,j are just constants. 

We can solve for u in the first line of the system to get u = (v -+ 1)/(j+1). Then substitute this into the two second order equations and solve. That's pretty much what the program is attempting to do.

Since the solve command in the program does actually find a solution, the eight conditions I placed on the solution should all hold, yet none of them do, so something I am doing is wrong.

Sorry, I think I must have changed the value of B and not ran it again when I posted. I apologize. Anyway, I get the same values as you guys now, but it doesn't change the fact that the values should be 0 if the conditions used to find the constants are true. 

The condition 

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

suggests that

vnegc(xi)-v1posc(xi) = 0, but from the code we obtain -3.221177116 (when xi = 0.2)

while the condition

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

suggests that 

vlnegc'(0) = 0, but we obtain -0.4814191895e-1 if we choose xi = 0.2 (and much larger if xi is larger)

Sorry again for the trouble. Any explanation of this issue would be greatly appreciated.

I actually can't upgarde at the current time. I will just solve it in two parts by doing the homogeneous first. Thanks for the help.