Since you want to avoid differential equations, I will instead consider the so-called Phasor domain, which is actually nothing but the Fourier transform of the original signals.
In phasor domain, we will be basically considering complex values: Complex voltages, complex resistances (which are denoted by $Z$ and called impedance): This is simply the mathematical convenience, and we will always get back the physical real values at the end.
For the time being, forget about the active parts of the circuit and focus on the passive elements, that is forget the switch and the DC source. Now let us check the behavior of remaining elements if we were to give the complex voltage $V(t)=e^{i\omega t}$:
For capacitor, we have the equation $C\frac{dV}{dt}=I$, hence we get $$I=i\omega C e^{i\omega t}$$ meaning the impedance is $$Z=\frac{V}{I}=\frac{1}{i\omega C}$$
For inductor, we have the equation $L\frac{dI}{dt}=V$, hence we get $$ i\omega L I=e^{i\omega t}$$ meaning the impedance is $$Z=\frac{V}{I}=i\omega L$$
Since the capacitor and the inductor are in series, the total impedance is their sum (just like the case for resistances), hence the impedance of the passive part of the circuit is $$Z=i\omega L+\frac{1}{i\omega C}=\frac{1-\omega^2 LC}{i\omega C}$$
This means that the voltage across the capacitor for a sinusoidal AC signal would be
$$V_\text{capacitor}=V_\text{input}\frac{Z_\text{capacitor}}{Z_\text{total}}=V_\text{input}\frac{\frac{1}{i\omega C}}{\frac{1-\omega^2 LC}{i\omega C}}$$
hence the ratio of input voltage to output voltage is
$$\frac{V_\text{output}}{V_\text{input}}=\frac{1}{1-\omega^2 LC}$$
where we took the output voltage to be the voltage of the capacitor.
One can immediately notice that the ratio would go to infinity if we were to excite the system with an AC source of the frequency $$\omega=\frac{1}{\sqrt{LC}}$$ This is called the resonance frequency. At this frequency, the system builds up each incoming power cycle and increases its output indefinitely. In real world, the existence of resistors in the systems and breakdown of equipment prevents this divergence to infinity.
If we now go back to our question, we now simply need to know how we can write our excitation in terms of pure sinusoidal waves of constant frequencies. If we know this, then we can examine the output of the system because we are dealing with a linear system: Effect of sum of sinusoidal signals is equal to sum of effects of individual signals.
Our input signal is simply step function, or so called Heaviside-theta function, which is zero before $t=0$ and constant after $t=0$ where we choose $t=0$ as the time the switch is turned on. One now needs to decompose Heaviside-theta function into pure sinusoidals: I'll skip the calculation here; it is basically Fourier transformation of this function. The answer is $$V_\text{input}(f)=E\left(\frac{i}{\sqrt{2\pi}f}+\sqrt{\frac{\pi}{2}}\delta(f)\right)$$
meaning that we are creating an input signal which decomposes of $V_\text{input}(f)$ for all frequencies $f$. We can now find the output signal for individual $f$ values by using the ratio above:
$$V_\text{output}(f)=\frac{E}{1-4 \pi^2 f^2 LC}\left(\frac{i}{\sqrt{2\pi}f}+\sqrt{\frac{\pi}{2}}\delta(f)\right)$$
where we used the fact that natural frequency $\omega$ is $2\pi f$.
We can now get our result taking the inverse Fourier transform, effectively summing over the contribution of all pure sinusoidals with different $\omega$ values: Again, this is due to the fact that the system is linear. I will skip the calculations, but you can easily use Mathematica to do this for you as can be seen at the attached picture. The result is
$$V_\text{output}(t)=\left\{\begin{aligned}E\left(1-\frac{1}{2} \cos \left(\frac{t}{2 \pi\sqrt{C L}}\right)\right)\quad t>0\\0\quad t<0\end{aligned}\right.$$
which is mathematically true but with the wrong boundary condition. I comment about this further down. Nonetheless, it shows us the relevant behavior.
You can see in the picture that the voltage simply keeps oscillating. The graphs reflect three facts:
- Energy keeps shifting from being stored in capacitor and being stored in the inductor.
- As there is no resistance, there is no damping effect, so the oscillations do not die!
- The output form is only dependent on the characteristics of the passive circuit; namely the inductance of inductor and the capacitance of capacitor. The input voltage only scales the output form. This is a generic concept: In engineering, this characteristic is called the impulse response whereas it is known to be the Green's function in physics, though impulse response is actually a specific case of the more general Green's functions.
EDITS (Answers to comments):
- There was a typo earlier where I wrote $\sqrt[4]{LC}$ instead of $\sqrt{LC}$. I corrected it above, everything is now dimensionally consistent.
The result respects the boundary value at $t=0$; however, this boundary value is ill-defined. We know that the input voltage is $E$ for $t>0$ and $0$ for $t<0$; however, we do not know what it is exactly at $t=0$ since turning on a switch is a discontinuous jump. The mathematically usual choice is to take it to be the half value, hence $E/2$ at $t=0$, this is indeed in the definition of Heaviside theta in Mathematica as well, and this is why we get $E/2$ for $t=0$ in our output! I was completely incorrect here, I realized it thanks to @Alfred
If one wants to define the behavior of the switch differently, one can repeat the calculations accordingly (actually using the differential equations would be a lot easier). The output has nothing to do with switch but with the boundary condition chose. However, with the chosen boundary condition, in the end, the static response $E$ in the output would remain the same whereas the jump response $-\frac{E}{2}\cos(\frac{t}{2\pi\sqrt{L C}})$ would scale to $-E\cos(\frac{t}{2\pi\sqrt{L C}})$, hence we would have
$$
V_\text{output}=V_\text{capacitor}=E\left(1-\cos(\frac{t}{2\pi\sqrt{L C}})\right)
$$
hence
$$I=C\frac{dV_\text{capacitor}}{dt}=\frac{E}{2\pi}\sqrt{\frac{C}{L}}\sin\left(\frac{t}{2\pi\sqrt{L C}}\right)$$

Further Comments
IMO, the simplest approach to the question is to use the differential equations. However, since OP asked us to avoid this, I tried to explain it by using the facts that electrical systems are linear systems and that eingenfunctions of linear systems are exponential functions; that is pure sinusoidals. I believe this is the conceptually closest approximation to using differential equations, though not equally rigorous. I will address some issues in this part regarding the approach above.
First of all, the inverse fourier transform of our output,
$$V_\text{output}(f)=\frac{E}{1-4 \pi^2 f^2 LC}\left(\frac{i}{\sqrt{2\pi}f}+\sqrt{\frac{\pi}{2}}\delta(f)\right)$$ is actually
$$V_\text{output}=\text{sgn}(t) \sin ^2\left(\frac{t}{4 \pi \sqrt{C L}}\right)+\frac{1}{2}$$
which has the behavior as the plot below:

There are a few commands in order:
- The plot naively looks non-causal, as there is nonzero signal even before the switch is turned on. However that signal is actually the resonant signal we mentioned above. As we discussed there, the effective resistance (impedance) for this signal is zero for the circuit, so that signal can never disappear in a generic solution: It has nothing to do with the switch being turned on or off! The effect of switch is purely shifting the output, in a continuous manner, which is causal.
- In differential equations language, the resonant signal is the homogeneous solution which is present in the system unless killed off by a boundary/initial condition. The effect of the switch is the particular solution, which is simply the shift of the output signal in a causal manner.
- In the EDITS part above, I got confused and incorrectly stated that the reason voltage is $1/2$ at $t=0$ is because of the convention of Heaviside-theta. This is not true. I believe it has to do with the fact that signals are assumed to converge to zero at minus infinity in Fourier Transform, hence the mean value of the sine should be 0 for $t<0$, forcing signal to be $1/2$ at $t=0$ to remain continuous with the $t>0$ part. If we instead impose the condition that it should be $0$ at $t=0$, we get the correct result for $t>0$, and irrelevant/unphysical signal for $t<0$.
- A way better method to do this calculation is to use Laplace domain instead of Fourier domain, where we consider Laplace transformation from $0$ to $\inf$. This would allow us to implement boundary condition directly in the Laplace domain and we would not get the unphysical part for $t<0$. One can also do that; I avoided it as I am not aiming for a rigorous or efficient method but a conceptually easy one, and Phasor space is easier to understand IMO.
To justify the explanations above, let us set up the differential equations:
\begin{align}
I_c=&I_l \Rightarrow\\
\frac{d I_c}{d t}=&\frac{d I_l}{d t}+\frac{a}{L} \Rightarrow\\
C\frac{d^2 V_c}{d t^2}=&\frac{V_l}{L}+\frac{a}{L} \Rightarrow\\
LC\frac{d^2 V_c}{d t^2}=&V_i-V_c+a \Rightarrow\\
\left(LC\frac{d^2}{d t^2}+1\right)V_c=&V_i+a
\end{align}
where $V_i$ is the input voltage.
The homogeneous solution is that of $$\left(LC\frac{d^2}{d t^2}+1\right)V_c=0$$ which is simply $$V_c^\text{homogeneous}=b\cos\left(\frac{t+c}{\sqrt{LC}}\right)$$
This is exactly the same signal we observe for $t<0$ in the plot above.
The correct procedure is actually to use Green's function with the input signal with appropriate boundary conditions: This is what EE engineers would call impulse response. Laplace domain is ideally suited for this job for linear systems. My explanation above on the other hand corresponds to finding particular solution for $V_i=E$ instead, combining it with homogeneous solution, imposing $V_E(t=0)=0$ and discarding $t<0$ part.