# Solving Equations with interpolation functions

I'm trying to solve a system of two equations, both of which use the same interpolating function. If the background would help, I'm trying to find the depletion width of a p-n junction after I've evaporated a dopant (donor) onto it and I've let it diffuse for a while. I'm trying a simple case first, where the semiconductor has a constant acceptor concentration of .8 atoms/cm^-3 and the donor concentration is some interpolation function that I got by solving the diffusion equation using the Crank-Nicolson method. I find the junction (where the two concentrations are equal) and then I solve the two simultaneous equations: 1.) Derivative of potential is zero at both boundaries 2.) Continuity of the scalar potential @ the junction.

Here's my code (just signed up for Mathematica SE, please bare with me):

Solve[{0.8*lp == (solution[ln, 5])*ln, (-solution[0, 5])*(ln)^2 == 0.8*(lp)^2}, {lp, ln}]


where 'solution' is the interpolation function (the donor concentration in the semicondcutor). However, when I try and solve this system of equations it gives me the error: "Solve::inex: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help."

I'm not sure if this is something that Mathematica can even do, but any help would be greatly appreciated.

Edit: Here's my code for how I got the interpolation function, 'solution

M = 100; Δt = 0.05; Δx = 0.1; junc = 3.16294;
ρ0[x_] := N[Erfc[x]];
ρ[0, n_] := ρ[1, n];
ρ[M, n_] := ρ[M - 1, n];
sol[0] = Table[ρ[j, 0] -> ρ0[j*Δx], {j, 1, M - 1}];
sol[n_] := sol[n] = Module[{vars, eqns}, vars = Table[ρ[j, n], {j, 1, M     - 1}];
eqns = Table[ρ[j, n] - ρ[j, n - 1] == (Δt*(ρ[j + 1, n] - 2*ρ[j, n] + ρ[j - 1, n] + ρ[j + 1, n - 1] - 2*ρ[j, n - 1] + ρ[j - 1, n - 1]))/(2*Δx^2),
{j, 1, M - 1}] /. sol[n - 1]; Solve[eqns, vars][[1]]]
solution = Interpolation[Flatten[Table[Table[{Δx*j - junc, n*Δt, ρ[j, n]}, {j, 0, M}] /. sol[n], {n, 0, 100}], 1]]

• Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Nov 16 '16 at 6:55
• – user9660 Nov 16 '16 at 6:59
• We can't run your code, because you haven't provided solution and junc. – Feyre Nov 16 '16 at 10:31
• Just edited with code for 'solution' interpolation function. – hijasonno Nov 16 '16 at 18:21
• Try using FindRoot. Using the most updated information I'm aware of the solution to the system is $\left(ln,lp\right) = \left(0,0\right)$. Try FindRoot[{0.8*lp == (solution[ln, 5])*ln, (-solution[0, 5])*(ln)^2 == 0.8*(lp)^2}, {{lp, 1}, {ln, 1}}] – Marchi Nov 16 '16 at 20:34