Predicting potential outcomes (POs) of interventions from observational data is crucial for decision-making in medicine, but the task is also challenging due to the fundamental problem of causal inference. Existing methods are commonly designed for conditional average treatment effect (CATE) estimation where POs are only auxiliary quantities; yet, in finite samples, such methods tend to underperform at predicting POs. In this paper, we propose a novel causal diffusion model for predicting POs. Specifically, we leverage a tailored conditional denoising diffusion model to learn complex distributions of POs, where we account for the causal structure of the task (i.e., adjust for the covariate shift across treated vs. non-treated individuals) through our orthogonal diffusion loss designed for the problem with finite samples. A particular strength of our method is that it is highly flexible: it can handle a wider range of settings from causal inference, such as (i) binary and continuous treatments, (ii) single and multiple outcomes, and (iii) different causal quantities (e.g., POs but also CATE). Across a wide range of experiments, we show that our method achieves state-of-the-art performance.