kill(all)$

Maxima : reaction  Formula  , 'version 1.0' ,part: '3a'
peter.vlasschaert@gmail.com, 30/10/2019

load(lrats)$

formula(v):=block([],
kaa:makelist(v[i][1],i,1,length(v)),
kbb:makelist(v[i][2],i,1,length(v)),
kbb1:makelist(if  v[i][2] =1 
    then  v[i][2] : "" else  v[i][2] : v[i][2],i,1,length(v)),
kcc1:maplist(concat,kaa,kbb1),
mpp:concat(kcc1[1],kcc1[2]),
for i:3 thru length(v) do(
mpp: concat(mpp,kcc1[i])),
    return(mpp)
);         

elements reaction
-------------------------

p:[Na,O,H,S];

products 
=======

1e) NaOH    =>  a

p1a:[[Na,1],[O,1],[H,1]];

p1aa:formula(p1a);

2e) H2SO4  => b

p2a:[[H,2],[S,1],[O,4]];

p2aa:formula(p2a);

reaction products
=============

3e)  H2O  => c

p1b:[[H,2],[O,1]];

p1bb:formula(p1b);

4e)  Na2SO4 => d

p2b:[[Na,2],[S,1],[O,4]];

p2bb:formula(p2b);

reaction

pp:a*p1aa + b*p2aa  = c*p1bb + d*p2bb;

How to find : a,b,c,d
================

rem :find number of element on each side

a1)   Na  : both side
----------------------------

products :
-------------

a1:[]$
for i:1 thru length(p1a)  do(
if p[1] =p1a[i][1]  then
a1:endcons(a*p1a[i][2],a1))$
a1;

reaction products :
-------------------------

aa1:[]$
for i:1 thru length(p2b)  do(
if p[1] =p2b[i][1]  then
aa1:endcons(d*p2b[i][2],aa1))$
aa1;

equation:
-------------

aa2:map("=",a1,aa1);

Na ,equation:
------------------

pp;

aa3:[p[1],aa2[1]];

a     : products               ,' Na'
2*d : reaction products ,'Na'

other , Na,O,H,S
=============

xx3:[p[1],aa2[1]];

bb3:[p[2],a+4*b = c+4*d];
cc3:[p[3],a+2*b =2* c];
dd3:[p[4],b=d];

find : a,b,c,d
==========

p1:[a=2*d,a+4*b = c+4*d,a+2*b =2*c,b = d];
p2:solve(p1,[a,b,c,d]);

p3:subst([%r1=1],p2);

p3a:flatten(p3);

reaction :coef  
==========

p4:subst(p3a,pp);