

function plot_normal_domain(a,b,psi,phi)

t=linspace(a,b,1000); t=t';
psi_t=feval(psi,t);
phi_t=feval(phi,t);

pts_bound_low=[t psi_t];
pts_bound_up=[flipud(t) flipud(phi_t)];
pts=[pts_bound_low; pts_bound_up; pts_bound_low(1,:)];

fill(pts(:,1),pts(:,2),'c');
hold on;