MMa is still in its infancy as far as solving pde's analytically. Solving this problem by separation of variables is fairly straightforward, although somewhat messy.
Clear["Global`*"]
I will simplify your psi1
phi1[x_] = a + b x + c x^2;
The pde
eqn2 = D[phi2[x, y], {x, 2}] + D[phi2[x, y], {y, 2}] == 0
The boundary conditions
bc3 = phi2[x, 0] == vs - phi1[x]; bc4 = phi2[x, l] == vs + vd - phi1[x]; bc5 = phi2[tch + d, y] == 0; bc6 = phi2[-d, y] == 0;
Assume
phi2[x_, y_] = X[x] Y[y] eqn2/phi2[x, y] // Expand (* X''[x]/X[x] + Y''[y]/Y[y] == 0 *)
If we move the Y part to the rhs we see that the X part is the negative of the Y part. Each side must be equal to a constant. Let the X part be sinusoidal. So we have.
xeq = X''[x]/X[x] == -alpha^2 yeq = Y''[y]/Y[y] == alpha^2
DSolve the xeq for X[x]
DSolve[xeq, X[x], x] // Flatten; X1 = X[x] /. % /. {C[1] -> c1, C[2] -> c2};
And the yeq
DSolve[yeq, Y[y], y] // Flatten; Y1 = Y[y] /. % /. {C[1] -> c3, C[2] -> c4};
Convert the above to trig
ExpToTrig[Y1] // Collect[#, {Sinh[alpha y], Cosh[alpha y]}] & (* (c3 - c4) Sinh[alpha y] + (c3 + c4) Cosh[alpha y] *) Y1 = % /. {c3 - c4 -> c3, c3 + c4 -> c4} (* c3 Sinh[alpha y] + c4 Cosh[alpha y] *)
We have the pieces. Put them together as
phi2[x_, y_] = X1 Y1;
To start, look at bc5
bc5 (* (c3 Sinh[alpha y] + c4 Cosh[alpha y]) (c1 Cos[alpha (d + tch)] + c2 Sin[alpha (d + tch)]) == 0 *) Solve[bc5, c2] (* {{c2 -> -c1 Cot[alpha (d + tch)]}} *) c2 = c2 /. %[[1]];
Now bc6
bc6 (* (c3 Sinh[alpha y] + c4 Cosh[alpha y]) (c1 Sin[alpha d] Cot[alpha (d + tch)] + c1 Cos[alpha d]) == 0 *) bc6[[1, 1]] // Collect[#, c1] & (* c1 (Sin[alpha d] Cot[alpha (d + tch)] + Cos[alpha d])*) %/c1; alphaeq = % == 0 (* Sin[alpha d] Cot[alpha (d + tch)] + Cos[alpha d] == 0 *)
The above determines values of alpha. It is a transcendental equation solved numerically, once constants of the problem are assigned. There are an infinite number of then and the solution becomes an infinite sum over the values of alphas. Ignore the sum for now.
phi2[x, y] // Simplify (* c1 Csc[alpha (d + tch)] Sin[ alpha (d + tch - x)] (c3 Sinh[alpha y] + c4 Cosh[alpha y]) *)
Combine constants. c1 is redundant
c1 = 1 phi2[x_, y_] = phi2[x, y] // Simplify
Now
bc3 (* c4 Csc[alpha (d + tch)] Sin[alpha (d + tch - x)] == -a - b x - c x^2 + vs *)
Use orthogonalities of the alphas to solve for c4.
c4*Csc[alpha*(d + tch)]*Integrate[Sin[alpha*(d + tch - x)]^2, {x, -d, tch + d}] == Integrate[(-a - b*x - c*x^2 + vs)*Sin[alpha*(d + tch - x)]^2, {x, -d, tch + d}]//Simplify; c4rule = Solve[%, c4] // Simplify;
Keep thing simple. Don't assign it yet.
bc4 (* Csc[alpha (d + tch)] Sin[ alpha (d + tch - x)] (c3 Sinh[alpha l] + c4 Cosh[alpha l]) == -a - b x - c x^2 + vd + vs *)
Now orthogonality to solve for c3.
Csc[alpha*(d + tch)]*(c3*Sinh[alpha*l] + c4*Cosh[alpha*l])* Integrate[Sin[alpha*(d + tch - x)]^2, {x, -d, tch + d}] == Integrate[(-a - b*x - c*x^2 + vd + vs)*Sin[alpha*(d + tch - x)]^2, {x, -d, tch + d}]//Simplify; c3rule = Solve[%, c3] // Simplify;
Now assign c3, c4
c4 = c4 /. c4rule[[1]] // Simplify; c3 = c3 /. c3rule[[1]] // Simplify;
We have solved for all the constants satisfying the boundary conditions which gives us phi2[x,y] in terms of alpha. The total solution is an infinite sum over all the alphas, determined by numerically solving the equation for alpha once numerical values for all your constants are specified. I have not printed the output for the later values of these calculations because they are lengthy.
