clear,clc,close all;
global FEM_M FEM_K FEM_A FEM_Q FEM_G FEM_F;
model = createpde(1);
R1 = [3;4;-1;1;1;-1;-1;-1;1;1];
g = decsg(R1);
heatmodel_geom = geometryFromEdges(model,g);
msh = generateMesh(model,'GeometricOrder','linear');
[P,E,T] = meshToPet( msh );
specifyCoefficients(model,'m',0,'d',1,'c',1,'a',0,'f',1);
applyBoundaryCondition ( model,'neumann','Edge',1:4 ,'g' ,-2,'q' ,3);
model_FEM_matrices = assembleFEMatrices(model);
FEM_M = model_FEM_matrices.M;
FEM_K = model_FEM_matrices.K;
FEM_A = model_FEM_matrices.A;
FEM_Q = model_FEM_matrices.Q;
FEM_G = model_FEM_matrices.G;
FEM_F = model_FEM_matrices.F;
ntime = 20;
startTime = 0;
endTime = 0.1;
tlist = linspace(startTime,endTime,ntime);
odesolve_tol = 1e-6;
u0 = IC_general(P(1,:), P(2,:));
options = odeset('Mass',FEM_M ,'AbsTol',odesolve_tol ,'RelTol',odesolve_tol ,'Stats','on');
disp('ode45');
[TOUT , YOUT] = ode45(@odefun_semidiscretize_pde,tlist,u0',options );
pdeplot(model ,'XYData',YOUT(end ,:),'ZData',YOUT(end ,:),'Contour','on','ColorMap','jet');
function f = IC_general(x,y)
nr = length(x);
f = zeros(1, nr);
f(1 ,:) = 1;
end
function Yout = odefun_semidiscretize_pde(t,Y)
global FEM_M FEM_K FEM_A FEM_Q FEM_G FEM_F
Yout = -(FEM_K*Y+ FEM_A*Y+ FEM_Q*Y)+ FEM_G + FEM_F ;
end