@tomleslie can we do this without initial/boundary conditions? Since they are not provided.

@tomleslie How to do next step? I want to integrate (4) for one more times so that we have first-order diff eq and then solve it by using 'dsolve'. 

@Christian Wolinski I want to extarct Sol1: 'Sol1={c = c0, mu = 1/sqrt(2*c0^2 - 2), a[-1] = a[-1], a[0] = 0, a[1] = a[1], b[-1] = a[-1], b[0] = 0, b[1] = -a[1]}'

sols_tw.mw

@Christian Wolinski Why we need to put 'map' and 'identity' in 'Sols'? Secondly, how can I extract 'c, mu, a[-1], a[0], a[1], b[-1], b[0], b[1]' from 'Sols'? By writing 'Sol1 := op(11, op(7, Sols));
Sol2 := op(15, op(7, Sols));' in the worksheet, I found an error.