function pdemodel
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',4);
pdetool('snapon','on');
set(ax,'DataAspectRatio',[1 1 1]);
set(ax,'PlotBoxAspectRatio',[1.5 1 1]);
set(ax,'XLim',[-1 2]);
set(ax,'YLim',[-0.5 1.5]);
set(ax,'XTickMode','auto');
set(ax,'YTickMode','auto');
pdetool('gridon','on');

% Geometry description:
pderect([0 1 0 1],'SQ1');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','SQ1')

% Boundary conditions:
pdetool('changemode',0)
pdesetbd(4,...
'dir',...
2,...
str2mat('1','0','0','1'),...
str2mat('0','0'))
pdesetbd(3,...
'dir',...
2,...
str2mat('1','0','0','1'),...
str2mat('0','0'))
pdesetbd(2,...
'dir',...
2,...
str2mat('1','0','0','1'),...
str2mat('0','0'))
pdesetbd(1,...
'dir',...
2,...
str2mat('1','0','0','1'),...
str2mat('0','0'))

% Mesh generation:
setuprop(pde_fig,'trisize',Inf);
setuprop(pde_fig,'Hgrad',1.3);
setuprop(pde_fig,'refinemethod','regular');
pdetool('initmesh')
pdetool('refine')
pdetool('refine')
pdetool('refine')
pdetool('refine')

% PDE coefficients:
pdeseteq(1,...
str2mat('2*((260)./(2*(1+(0.3))))+(2*((260)./(2*(1+(0.3)))).*(0.3)./(1-2*(0.3)))','0','(260)./(2*(1+(0.3)))','0','(260)./(2*(1+(0.3)))','2*((260)./(2*(1+(0.3)))).*(0.3)./(1-2*(0.3))','0','(260)./(2*(1+(0.3)))','0','2*((260)./(2*(1+(0.3))))+(2*((260)./(2*(1+(0.3)))).*(0.3)./(1-2*(0.3)))'),...
str2mat('0.0','0.0','0.0','0.0'),...
str2mat('40000*(x-y).*(6*x.^3.*y - 3*x.^3 + 4*x.^2 -6*x.^2.*y.^2 -3*x.^2.*y+ 6*x.*y.^3-2*x.*y -3*x.*y.^2 -1 +4*y.^2 - 3*y.^3)','(-40000*(x-1+y).*(6*x.^3.*y - 3*x.^3 -15*x.^2.*y +5*x.^2+ 6*x.^2.*y.^2 -x +10*x.*y -15*x.*y.^2 + 6*x.*y.^3 + 5*y.^2 -y -3*y.^3))'),...
str2mat('1.0','0','0','1.0'),...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setuprop(pde_fig,'currparam',...
['260                                                                                                                             ';...
'0.3                                                                                                                             ';...
'40000*(x-y).*(6*x.^3.*y - 3*x.^3 + 4*x.^2 -6*x.^2.*y.^2 -3*x.^2.*y+ 6*x.*y.^3-2*x.*y -3*x.*y.^2 -1 +4*y.^2 - 3*y.^3)            ';...
'(-40000*(x-1+y).*(6*x.^3.*y - 3*x.^3 -15*x.^2.*y +5*x.^2+ 6*x.^2.*y.^2 -x +10*x.*y -15*x.*y.^2 + 6*x.*y.^3 + 5*y.^2 -y -3*y.^3))';...
'1.0                                                                                                                             '])

% Solve parameters:
setuprop(pde_fig,'solveparam',...
str2mat('0','1000','10','pdeadworst',...
'0.5','longest','0','1E-4','','fixed','Inf'))

% Plotflags and user data strings:
setuprop(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]);
setuprop(pde_fig,'colstring','');
setuprop(pde_fig,'arrowstring','');
setuprop(pde_fig,'deformstring','');
setuprop(pde_fig,'heightstring','');

% Solve PDE:
pdetool('solve')