(* Ejemplo de módulos *)
fa[x_,y_]:=Module[{a,sol},
sol=Solve[{x a+y b==1,x a-y b==2},{a,b}];
a/.sol[[1]]
]
fb[x_,y_]:=Module[{a,sol},
sol=Solve[{x a+y b==1,x a-y b==2},{a,b}];
b/.sol[[1]]
]
f[x_,y_]:=Module[{a,b,sol},
sol=Solve[{x a+y b==1,x a-y b==2},{a,b}];
{a/.sol[[1]],b/.sol[[1]]}
]

f[5,5]

sol

(* Limpiamos las variables del ejemplo anterior *)
Clear[x,y,f]

punto[punto1_,punto2_,coef_]:=Module[
{x1,y1,u1,p1,q1,x2,y2,u2,p2,q2,x3,y3,u3,p3,q3,x,y,u,p,q,cond1,cond2,f1,g2,a1,a2,b1,b2,c1,c2,e1,e2,f,g,sol,a,b,c,e},
(* punto1->{x1,y1,u1,p1,q1}
punto2->{x2,y2,u2,p2,q2} *)
a=coef[[1]];
b=coef[[2]];
c=coef[[3]];
e=coef[[4]];
x1=punto1[[1]];
y1=punto1[[2]];
u1=punto1[[3]];
p1=punto1[[4]];
q1=punto1[[5]];
x2=punto2[[1]];
y2=punto2[[2]];
u2=punto2[[3]];
p2=punto2[[4]];
q2=punto2[[5]];
cond1={x->x1,y->y1,u->u1,p->p1,q->q1};
cond2={x->x2,y->y2,u->u2,p->p2,q->q2};
f=((b+(Sqrt[(b^2)-a c]))/a);
g=((b-(Sqrt[(b^2)-a c]))/a);
f1=f/.cond1;
g2=g/.cond2;
a1=a/.cond1;
a2=a/.cond2;
c1=c/.cond1;
c2=c/.cond2;
e1=e/.cond1;
e2=e/.cond2;
sol=NSolve[{y-y1==(x-x1)f1,y-y2==(x-x2)g2},{x,y}];
x3=x/.sol[[1]];
y3=y/.sol[[1]];
sol=NSolve[{e1(x3-x1)f1-a1(p-p1)f1-c1(q-q1)==0,e2(x3-x2)g2-a2(p-p2)g2-c2(q-q2)==0},{p,q}];
q3=q/.sol[[1]];
p3=p/.sol[[1]];
u3=(1/2)(u1+(1/2)(x3-x1)(p1+p3)+(1/2)(y3-y1)(q1+q3)+u2+(1/2)(x3-x2)(p2+p3)+(1/2)(y3-y2)(q2+q3));
{x3,y3,u3,p3,q3}
]

(* Ejemplo de un par de puntos *)
punto1={0.0,0.0,0.0,0.0,0.0};
punto2={0.1,0.0,0.0,0.0,0.9708055193627332};
coef={1,0,-1,0};
punto[punto1,punto2,coef]

(*Ahora hacemos una línea completa *)
(* Primero calculo que ocurre en las fronteras *)
ue = -(1/2)Cos[π(x+y)]+(1/2)Cos[π(x-y)];
pe=D[ue,{x,1}];
qe=D[ue,{y,1}];
n=100;
puntos=Table[{1.0i/n,0.0,ue/.{x->1.0i/n,y->0.0},pe/.{x->1.0i/n,y->0.0},qe/.{x->1.0i/n,y->0.0}},{i,0,n}];

puntos
MatrixForm[puntos]

coef={1,0,-1,0};
punto[puntos[[1]],puntos[[2]],coef]

nd=Dimensions[puntos][[1]]

coef={1,0,-1,0};
plot={};
Do[
dato=punto[puntos[[i]],puntos[[i+1]],coef];
x3=dato[[1]];
u3=dato[[3]];
AppendTo[plot,{x3,u3}];
,{i,1,nd-1}]

plot
MatrixForm[plot]

ListPlot[plot]

Plot[ue/.{y->0.05},{x,0,1}]

Show[{Plot[ue/.{y->0.005},{x,0,1}],ListPlot[plot]}]