How can I model this linear power system in MATLAB?
To solve this problem you should look at the ODE solvers in MATLAB.  

Note that your function does not describe what the response of the system is to a variation in input power Pm.

Below I have written up an implementation, based on the example in http://www.its.caltech.edu/~ae121/ae121/ode45_Ref2.pdf .  In my example I have added a term Pm to the differential equations, so you can see how a disturbance might affect the dynamics.

Unfortunately I don't have MATLAB installed so I can't check whether there are any errors in this, but the idea should be clear.

function main 
% a simple example to solve ODE's
% Uses ODE45 to solve 
%    dx_dt(1) = 1*x(1)+2*x(2)+3*x(3)+4*x(4)
%    dx_dt(2) = 5*x(1)+6*x(2)+7*x(3)+8*x(4)+Pm
%    dx_dt(3) = 9*x(1)+1*x(2)+2*x(3)+3*x(4)
%    dx_dt(4) = 4*x(1)+5*x(2)+6*x(3)+7*x(4)
%set an error
options=odeset('RelTol',1e-6);

%initial conditions
X0 = [0;0;0;0]; 
Pm0=1;

%before disturbance
tspan1 = [-1,0]; 
[t1,X1] = ode45('system',tspan1,X0,options,Pm0);

%during disturbance
tspan2 = [0,1];
[t2,X2] = ode45('system',tspan2,X1(end),options,1.3*Pm0);

%after disturbance
tspan3 = [1,3];
[t3,X3] = ode45('system',tspan3,X2(end),options,Pm0);

time=[t1,t2,t3];
X=[X1,X2,X3];

%plot the results
figure 
hold on
plot(time,X)
legend('x1','x2','x3','x4');ylabel('x');xlabel('t') 
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dx_dt]= system(t,x,Pm) 
%a function which returns a rate of change vector
A = [1,2,3,4;
     5,6,7,8;... 
     9,1,2,3;
     4,5,6,7]
P=[0;Pm;0;0]
dx_dt = A*x+P; 
return