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}
]
linea1[puntosin_,coef_]:=Module[{nd,puntosout,dot},
nd=Dimensions[puntosin][[1]];
puntosout={};
Do[
dot=punto[puntosin[[i]],puntosin[[i+1]],coef];
AppendTo[puntosout,dot];
,{i,1,nd-1}];
puntosout
]
linea2[puntosin_,froni_,frond_,coef_]:=Module[{nd,puntosout,dot},
nd=Dimensions[puntosin][[1]];
puntosout={froni};
Do[
dot=punto[puntosin[[i]],puntosin[[i+1]],coef];
AppendTo[puntosout,dot];
,{i,1,nd-1}];
AppendTo[puntosout,frond];
puntosout
]

(* Calculamos x,y,u,p y q sobre la frontera *)
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}];

(* Este es el programa que resuelve la ecuación diferencial *)
nx=100;
ny=50;
plot={};
coef={1.0,0.0,-1.0,0.0};
Print["j=",Dynamic[j],"/",ny]
progj=0;
ProgressIndicator[Dynamic[progj]]
Do[
progj=1.0 j/ny;
puntos1=linea1[puntos,coef];
(* Calcular los valores en la frontera *)
npuntos1=Dimensions[puntos1][[1]];
puntoi=punto[puntos1[[1]],puntos1[[2]],coef];
puntod=punto[puntos1[[npuntos1-1]],puntos1[[npuntos1]],coef];
xi=0.0;
yi=puntoi[[2]];
ui=ue/.{x->xi,y->yi};
pi=pe/.{x->xi,y->yi};
qi=qe/.{x->xi,y->yi};
xd=1.0;
yd=puntod[[2]];
ud=ue/.{x->xd,y->yd};
pd=pe/.{x->xd,y->yd};
qd=qe/.{x->xd,y->yd};
froni={xi,yi,ui,pi,qi};
frond={xd,yd,ud,pd,qd};
puntos2=linea2[puntos1,froni,frond,coef];
puntos=puntos2;

nx=Dimensions[puntos][[1]];
(* Aquí generamos un archivo que contiene 
las coordenadas x,y de cada punto de la red y el
valor de u en cada uno de estos puntos*)
Do[
x3=puntos[[i,1]];
y3=puntos[[i,2]];
(* u3=puntos[[i,3]]; *)
p3=puntos[[i,4]];
AppendTo[plot,{x3,y3,p3}];
,{i,1,nx}]
,{j,1,ny}]

ListPlot3D[plot]

Plot3D[pe,{x,0,1},{y,0,0.45}]

Dimensions[plot]

out=plot[[3543]]

val=pe/.{x->out[[1]],y->out[[2]]}

100(val-out[[3]])/val

patito[a_,x1_,y1_]:=Module[{},
a/.{x->x1,y->y1}
]

patito[1,2,2]

patito[(x^2)+y,3,4]