Modeling chemical reactions using system of equations

There are 4 proteins that form a reaction network as follows: P1 can change into P2 with rate constant k12. This reaction is reversible with a constant k21. P2 can be changed into P3 and P4 with rate constants k23 and k24 respectively. These reactions are also reversible with rate constants k32 and k42. At t=0 there is only P1 present in the solution and t=0 P1=1

%Rate constants
k12 = 0.5;
k21 = 0.1;
k23 = 0.3;
k24 = 0.7;
k32 = 0.2;
k42 = 0.2;

delta = 0.01;%increment for x variable

I know that in order to solve this system I have to create a coefficient matrix, create a vector for x, create a dependent array, and then fill the rest of it using Euler's formula. I'm just not sure how to code everything.